Accurate, efficient and simple forces with Quantum Monte Carlo methods 
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Computation of ionic forces using quantum Monte Carlo methods has long been a challenge. We 
introduce a simple procedure, based on known properties of physical electronic densities, to make 
the variance of the Hellmann-Feynman estimator finite. We obtain very accurate geometries for the 
molecules H2, LiH, CH4, NH3, H2O and HF, with a Slater- Jastrow trial wave function. Harmonic 
frequencies for diatomics are also in good agreement with experiment. An antithetical sampling 
method is also discussed for additional reduction of the variance. 



The optimization of molecular geometries and crys- 
tal structures and ab initio molecular dynamics simula- 
tions are among the most significant achievements of sin- 
gle particle theories. These accomplishments were both 
possible thanks to the possibility of readily computing 
forces on the ions within the framework of the Born- 
Oppcnhcimer approximation. The approximate treat- 
ment of electron interactions typical of these approaches 
can, however, lead to quantitatively, and sometimes qual- 
itatively, wrong results. This fact, together with a favor- 
able scaling of the computational cost with respect to 
the number of particles, has spurred the development of 
stochastic techniques, i.e. quantum Monte Carlo (QMC) 
methods. Despite the higher accuracy achievable for 
many physical properties, the lack of an efficient estima- 
tor for forces has prevented, until recently 0,0,0! the use 
of QMC methods to predict even the simplest molecular 
geometry. The chief problem is to have a Monte Carlo 
(MC) estimator for the force with sufficiently small vari- 
ance. For example, in all-electron calculations, a straight- 
forward application of MC sampling of the Hellmann- 
Feynman estimator has infinite variance. This can be 
easily seen from the definition of the force. For a nu- 
cleus of charge Z at the origin, the force can be written, 
together with its variance, as a function of the charge 
density p(r) as 




F = Z drp(r 



F 



a = Zj 



dvp{v)^-F 2 . (1) 



Since the electronic density is finite at the origin, the 
variance integral diverges. 

In this paper, we propose a modified form for the 
force estimator which has finite variance. This estima- 
tor is then used to calculate forces and predict equilib- 
rium geometry and vibrational frequencies for a set of 
small molecules. Without loss of generality we will con- 
sider only the z-componcnt of the force on an atom at 
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FIG. 1: Force density along the z-direction for the H atom in 
LiH. The bond is along the z-axis, with a length of 3.316 Bohr. 
The continuous black curve is calculated from the Hartree- 
Fock orbitals. The dashed line is the estimate of f z using 
the bare estimator. Circles are obtained in an identical QMC 
simulation using the antithetic sampling technique outlined 
in the text. 



the origin. In a QMC calculation based in configura- 
tion space, the charge density is a sum of delta func- 
tions: p(r) oc J2 r > — r '), where the sum is over all 
N e electron positions and all MC samples. We consider- 
separately the electrons within a distance 1Z of the atom 
and those outside. The contribution to the force from 
charges outside, F®, can be calculated directly with the 
Hellmann-Feynman estimator in Eq. The contribu- 
tion from inside the sphere is responsible for the large 
variances in the direct estimator. It is convenient to in- 
troduce a "force density" defined as the force arising from 
electron charges at a distance r from the origin: 



f z (r)=Z / dn P (r,0 



I cos ( 



(2) 



Then the force is given as: 
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f z {r)dr. 



(3) 



The force density is a smooth function of r that tends to 
linearly as r approaches the origin. The force density for 
H in a LiH molecule computed with Hartree-Fock and two 
different QMC estimators is shown in Fig^ As expected 
the bare force estimator fluctuates wildly at small r. 

Because the force density is a smooth function, we can 
represent it in the interval (0, TV) with a polynomial 
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and determine the coefficients, a^, by minimizing 



x 2 = 



dr i 



f*(r)-fz(r) 



(4) 



(5) 



where r m is a weight factor used to balance contributions 
from different values of r. 

Since the relation between the force and the force den- 
sity is linear, and the relation between the fitting coeffi- 
cients and the electronic density is linear, we can directly 
write the force as averages over moments of the force den- 
sity. After some manipulations we arrive at: 



MC 



where the new estimator function is: 

M 

g(r) = e(K-r)J2c k r k+m . 



(6) 



(7) 



k=l 



The coefficients Cfc's are determined by c = S -"-h where 

(8) 



the Hilbcrt matrix S and the residual vector h are 

j^in+k+j+l TV +1 
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FIG. 2: Dependence of the VMC force on the expansion basis, 
for LiH with a bond length of 3.316 Bohr. The fitting radius 
1Z =0.6 Bohr. The definitions of the basis functions are in 
Eq.'s and @. The forces on H and Li are different because 
of the lack of full optimization of the VMC wave function (see 
text). 
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FIG. 3: Projection of the force in LiH using forward walking. 
The points at negative imaginary time give the VMC values. 
Values at are the mixed estimates of the DMC simulation. 



Note that for the bare estimator g(r) = 8(71 — r). Be- 
cause of the restriction on the basis, the variance of the 
new estimator is finite as long as m > —1/2. We have nu- 
merically found that the weighting factor m = 2, where 
each volume element is weighted equally, gives the lowest 
variance estimate of the force. 

To derive the estimator we have used the fact that 
f z (r) goes linearly at small r. 01 This is the crucial 
property that allows to filter out the s-wave component 
of the density responsible for the variance divergence. 
The original estimator is correct for any arbitrary charge 
density while the new filtered one uses physical proper- 
ties of the charge density to reduce the variance. The 
variance depends on the fitting radius 1Z and on the ba- 
sis set size M. As 1Z increases, the size of the basis must 
increase, which increases the variance. Charge densities 
corresponding to low energy states must be smooth and 
we typically find that only 2 or 3 basis functions are 
needed. The size of the basis can be reduced by using 
more appropriate basis sets. For example, in all calcula- 
tions reported below we used the expansion 



M 

fc=0 



(9) 



where /| D is the force density of a single determinant 
wave function, which can be readily computed from the 
orbitals. The improved basis allows a smaller polyno- 
mial set and a reduction of the variance. In Fig. [2] the 
dependence of the bias on the basis set type and size is 
shown for the case of a variational Monte Carlo (VMC) 
simulation on LiH at a bond length of 3.316 Bohr. 

The trial wave functions used in all cases were of 
the Slater-Jastrow form. The orbitals were obtained from 
a Hartree-Fock calculation using CRYSTAL98 0. The 
electron-electron and electron-proton Jastrow factors had 
the form of exp(ar/(l + br)), with a and b optimized 
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TABLE I: Equilibrium distances in A. Experimental, 
CCSD(T) and B3LYP values were taken from Ref@. The 
CCSD(T) and the B3LYP results were obtained using the cc- 
pVTZ basis set with the exception of LiH where the 6-311G* 
set was used. PBE results [3 were all obtained using the 
aug-cc-pVTZ basis set. 





QMC 


Exp. 


CCSD(T) 


B3LYP 


PBE 


H 2 


0.7419(4) 


0.741 


0.743 


0.743 


0.751 


LiH 


1.592(4) 


1.596 


1.618 


1.595 


1.606 


CH4 


1.091(1) 


1.094 


1.089 


1.088 


1.096 


NH 3 (N-H) 


1.009(2) 


1.012 


1.014 


1.014 


1.023 


NH 3 (H-H) 


1.624(2) 


1.624 


1.616 


1.624 


1.634 


H 2 (O-H) 


0.959(2) 


0.956 


0.959 


0.961 


0.971 


H 2 (H-H) 


1.519(3) 


1.517 


1.508 


1.520 


1.531 


HF 


0.919(1) 


0.918 


0.917 


0.923 


0.932 



by minimizing \Ei oc — (E)\ Q over points sampled from 
|* T | 2 - The time step in the diffusion Monte Carlo (DMC) 
simulations was chosen to give an acceptance ratio of 
98%, a value for which the time step bias on forces was 
within the statistical error bars. 

Since the exact density is needed for the Hellmann- 
Feynman theorem, forward walki ng N| or one of the varia- 
tional path integral algorithms|9lll0j] is needed in order to 
evaluate the force estimator. An example of the conver- 
gence of forward walking is shown in Fig. 01 The force as 
a function of the forward-walking projection time quickly 
reaches a plateau corresponding to the exact value. In 
this example, the variational forces are far from correct. 
This discrepancy results from the lack of full optimization 
of the trial wave function made of localized basis orbitals 
and atom centered Jastrow factors, and can be reduced 
somewhat by including Pulay terms Q ■ In DMC, forward 
walking eliminates the need for the Pulay corrections. 

The equilibrium geometries were computed by fitting 
the QMC forces in the proximity of the equilibrium ge- 
ometry to a polynomial with the appropriate symmetry. 
Fig. 0] shows the force in hydrogen fluoride in a 2% in- 
terval around the equilibrium geometry. The equilib- 
rium geometries are reported in Tabic ^ together with 
those given by CCSD(T), DFT using the B3LYP and 
the PBE functional, and experiments. The differences 
between QMC and experimental values are in all cases 
less than 0.4% and closer to the experiment than the 
other techniques. For diatomics it is easy to provide an 
estimate of the harmonic vibrational frequencies starting 
from the derivative of the force curve at equilibrium ge- 
ometry. The QMC frequencies, reported in Table |nl arc 
in good agreement with the experiment, with errors com- 
parable to that from CCSD(T) and DFT PBE or B3LYP. 
This suggests that forces computed within our approach 
are accurate also away from the equilibrium and could be 
used in molecular dynamics calculations or to optimize 
molecular geometries. 

The only source of systematic error in our calculations 



TABLE II: Harmonic frequencies in cm . Experimental, 
CCSD(T) and B3LYP values were taken from Ref@. The 
CCSD(T) and the B3LYP results were obtained using the cc- 
pVTZ basis set with the exception of LiH where the 6-311G* 
set was used. PBE results 114 were obtained using ad hoc 
gaussian basis sets. 





QMC 


Exp. 


CCSD(T) 


B3LYP 


PBE 


H 2 


4464(18) 


4410 


4420 


4401 


4323 


LiH 


1445(20) 


1369 


1414 


1405 


1380 


HF 


4032(266) 


4181 


4085 


4138 


4001 



that cannot be simply addressed is the fixed-node er- 
ror. In fixed-node DMC, the random walk is forbidden to 
cross the nodes of the trial wavefunction in order to pre- 
vent the loss of efficiency due to the fcrmion antisymme- 
try. If the nodes are accurate, so is the QMC energy and 
electronic density; hence the force. For incorrect nodes, 
the energy is an upper bound to the true energy, but such 
can not be said for the force. It is also not necessarily the 
case that the forces obtained from Eg. (1) are equal to 
the gradient of the fixed-node energy[ll],0,0|: this is 
only guaranteed in the limit of exact nodal surfaces. The 
high quality of the geometries and vibrational frequencies 
suggests that these errors, at least for the cases treated 
in this paper, are negligible. This is perhaps not surpris- 
ing, since the electronic density is a 1-electron property, 
while the nodal error is a many-body effect. 

We have also tested another method to further reduce 
the variance of the Hcllmann-Fcynman estimator. The 
filtered estimator performs well on the hydrogen atom 
but for heavier nuclei the error bar grows and seems to 
scale as Z 3 . In those cases the new method can poten- 
tially be very useful, with error bars scaling between Z 
and Z 2 . The method is based on the observation that, 
while electrons in the core cause large fluctuations in the 
force density, they contribute very little to it. A stan- 
dard approach to reduce the variance of a Monte Carlo 
estimate is the use of antithetic variates|l5|: a positive 
fluctuation is paired with a negative fluctuation. Sup- 
pose the random walk arrives at a multidimensional elec- 
tronic configuration R, with p (> 1) electrons inside a 
radius lZ av < 1Z of an atom located at the origin. We ob- 
tain an antithetic configuration R' by reflecting all p core 
electrons about the origin. Wc then estimate the force 
contribution due to the p electrons using both R and i?', 
assigning a weight factor of w(R') = \ip(R')/ip(R)\ 2 to 
R' . Their joint contribution to the estimator in Eq. © 

is Z- — W £ R ' g( r i) z i r ~ 3 where the sum runs over the p 
core electrons. Since w — > 1 as TZ av —> 0, fluctuations in 
the core are much reduced. 

Within VMC this scheme can be implemented exactly, 
leading to a dramatic reduction of the variance as can be 
noticed from Fig. ^ However this estimator is non-local 
and, in DMC, suffers from the same problems as non- 
local pscudopotentials, making an unbiased implcmcnta- 
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FIG. 4: DMC force in hydrogen fluoride. Left panel: evo- 
lution of the force over forward-walking time. Right panel: 
fully projected forces as a function of nuclear distance. R eq is 
the experimental equilibrium distance. 



tion not straightforward. We postpone further discussion 
of the antithetic method to a future article. 

Two other approaches have been introduced recently 
for the computation of forces in QMC. Filippi and Um- 
rigar have computed forces for diatomics by correlating 
random walks for interatomic separations a and a'. In 
DMC the difficulty associated with the nodal error and 
the branching factor was overcome by neglecting some 
types of correlation. The main drawback of a finite dif- 
ference method is the difficulty of calculating all the com- 
ponents of the force simultaneously; for a system of N 
atoms this method would require 3N separate force cal- 
culations. 

The other approach, introduced in Ref . Q , is closer to 
our method. It is based on a "zero-variance" version of 



the Hcllmann-Fcynman estimator and can be understood 
in the framework of this paper: one can prove that it 
corresponds to filtering out the s-wave component of the 
density leaving the force density unchanged. The semi- 
local character of the "zero- variance" estimator makes its 
DMC implementation trickier. To overcome this problem 
there have been attempts @, E3 to use correction terms 
similar in nature to the Pulay terms in single-particle ap- 
proaches. In practice, this scheme requires extensive op- 
timization and, although promising, it is unclear if it will 
be viable for more complicated cases. In addition, the 
value of the force is very sensitive to small errors |l6| in 
the charge density and the optimization within a stochas- 
tic technique is probably not sufficiently stable to elimi- 
nate these errors. 

In conclusion, we have developed a simple method for 
computing forces within quantum Monte Carlo and used 
it to find the equilibrium geometries for small polyatomic 
molecules. This has been the first time that a QMC 
technique is used to predict geometries of molecules be- 
yond diatomics. The only overhead in the calculation is 
the necessity of determining unbiased estimators, which 
requires the use of either forward-walking or reptation 
MC techniques. The new method leads to very accu- 
rate forces despite errors from the fixed-node approxima- 
tion and from its contribution to the energy derivatives. 
Extension of the method, including the antithetic esti- 
mator technique . to heavier atoms and to atoms with 
pseudopotentials |l8j is under investigation. 
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The force density f z (r) is proportional to the p z compo- 
nent of the density. A non-zero value as r — > would im- 
ply a discontinuity of p at the origin along the z-direction. 
The filtered estimator can be applied to atoms with non- 
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local pseudopotentials, where the density is replaced by be defined by summing over partial waves, 

a density matrix and a corresponding force density can 



